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Abstract 

A multifunctional interface method with 
capabilities for variable-fidelity modeling and multiple 
method analysis is presented. The methodology 
provides an effective capability by which domains with 
diverse idealizations can be modeled independently to 
exploit the advantages of one approach over another. 
The multifunctional method is used to couple 
independently discretized subdomains, and it is used to 
couple the finite element and the finite difference 
methods. The method is based on a weighted residual 
variational method and is presented for two- 
dimensional scalar-field problems. A verification test 
problem and a benchmark application are presented, 
and the computational implications are discussed. 

Introduction 

Accurate response characterization is needed early 
in the design and analysis cycle to investigate rapidly 
novel and revolutionary aerospace and ground vehicle 
design concepts. Mathematical modeling 

approximations are considered that range from simple 
handbook equations, empirically derived relations, 
spreadsheets, and design charts to complex continuous 
and discrete simulation models. Integrated 
computational methods with multiple capabilities and 
with diverse engineering applications provide the 
enabling technologies to assimilate variable-fidelity, 
multi-disciplinary data rapidly at an early stage in the 
design. These methods allow rapid design trade-offs 
between cost and performance. Based on the insight 
provided by these simulations, design uncertainties and 
risk assessment may be evaluated. Moreover, 
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methodology that allows different mathematical 
modeling approximations of the physical phenomena 
and among multiple engineering disciplines can be very 
advantageous. 

One class of methods for obtaining the desired 
modeling flexibility in the analysis of complex systems 
that has received recent attention is generally called 
interface methods. These methods, through the 
subdivision of a single domain into multiple 
subdomains, are used to couple independently 
discretized models or to couple different mathematical 
modeling methods or approximations, and exploit the 
advantages of one approach over another. Examples of 
methods that couple different independently discretized 
models include global/local approaches 1 * " 4 and finite 
element coupling methods 5 " 10 . Examples of methods 
that couple different modeling methods include finite 
element and boundary element methods 11 " 14 , finite 
element and Rayleigh-Ritz approximations 15 , finite 
element and finite difference methods 16 ' 17 , finite 
element and analytical solutions 18 , and finite element 
and equivalent plate solutions 19 . These interface 
methods facilitate the discretization of geometry by 
providing a capability to model independently the 
regions of interest, increasing the discretization fidelity 
or enhancing the mathematical approximation only in 
desired domains. A major feature of the methods is the 
accurate mapping of field quantities across the 
respective interface, whether the interface is one among 
diverse mathematical approximations or among diverse 
disciplines. Many of the interface methods are specific 
to a single spatial modeling approach ( i.e ., finite 
element method) or discipline, which limits their 
applicability to a wide range of applications in 
engineering science. 

This paper presents a multifunctional interface 
method with capabilities for variable-fidelity modeling 
and multiple method analysis. An interface method 
referred to in the literature as interface technology 9 ' 10 
was developed for coupling independently discretized 
finite element models and is well suited for extensions 
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to multiple method and multidisciplinary analysis. 
Therefore, the interface technology is the concept upon 
which the multifunctional capability presented in this 
work is built. The basic concept of the interface 
technology was discussed by Housner and Aminpour 9 
and subsequently, developed, implemented and 
validated by Aminpour et al. 10 , and Ransom et al. 2n . 
The method has been shown to be computationally 
efficient while preserving solution accuracy, however, 
its demonstration has been limited to solid mechanics 
problems. 

The purpose of this paper is to investigate 
multifunctional interface methods, as described herein 
that address the engineering design and analysis needs 
of multidisciplinary problems in engineering science. 
Fundamental relationships will be discussed among 
underlying engineering science and mechanics 
principles, computational methods, variable-fidelity 
models, and multiple methods using basic two- 
dimensional scalar-field problems from continuum 
mechanics. A multifunctional approach for scalar-field 
problems is developed, described, and demonstrated on 
problems in engineering science. A multifunctional 
formulation is presented along with a discussion of the 
spatial modeling and the computational implications. 
The multifunctional method is used to couple 
independently discretized subdomains and to couple the 
finite element and the finite difference methods. 
Numerical results for a verification test case and a 
representative scalar-field problem are presented. 

Overview of the Multifunctional Approach for 
Scalar-field Problems 


Methods of approximation such as Ritz, Galerkin, 
and other weighted residual methods are based on weak 
statements of the differential equations governing the 
system response. The scalar-field problem is a basic 
form of the governing differential equations; and thus, 
it lends itself to forming the mathematical foundation 
for the general methodology developed herein. 
Representative examples of scalar-field differential 
equations in two dimensions are considered, and the 
mathematical statement of the multifunctional approach 
is formulated. 

The general form of the differential equation 
describing a scalar-field problem for domain £2 (see 
Figure 1(a)) is given by the Poisson equation, which is 
of the form 


- V ■ {kSJu) = Q in £2 (1) 


subject to the natural boundary condition, 

k — + h(u ~u oa ) = q on T s , and essential boundary 
d« 


condition, u = it on Y p The normal 

, . . d u du du 

derivative, — = — n,. -I n„, and n. and n v are the 

dn dx x dy y y 

components of the outward normal vector, n, to the 
bounding surface, T, of domain, £2. In Eq. (1), the 
variables k and Q are known coefficients, and the 
primary variable or dependent variable is n, which is a 
function of the independent variables, x and y. In the 
natural boundary condition, the variables, h and are 
the convection coefficient, and the far-field value of the 

du 

primary variable, respectively. The terms, q, k — , and 

dx 

k — are the secondary variables that may be described 
dy 

on a portion of the boundary, T s . The primary 
variable, n, is specified on the boundary, T 7 ', and its 
prescription to the boundary value, u , constitutes the 
essential boundary condition. The complete boundary 
is defined as T = T p + T* . 



(a) Geometric representation 



(b) Boundary Definitions for subdomains 


Figure 1 . Geometric Representation and Boundary 
Definitions of Two-Dimensional Domain. 


Equations of the type of Eq. (1) arise in many 
fields of engineering science such as elasticity, heat 


2 


American Institute of Aeronautics and Astronautics 


transfer, fluid mechanics, and electrostatics. Reddy 21 
has tabulated several examples of engineering fields in 
which the analytical description of the physical process 
is described by Eq. (1). In this paper, the Poisson 
equation is applied to selected problems in the solid 
mechanics. 

Multifunctional Formulation 

The method of weighted residuals is used 
extensively in fluid mechanics; and thus, the potential 
problem is formulated from this perspective. The 
general weighted residual formulation for a single 
domain is presented in reference 22. The 
multifunctional formulation provides a capability to 
subdivide a domain into multiple, independently 
idealized subdomains, which allows the use of 
independently generated models, multiple levels of 
spatial discretization, and different mathematical 
approximations among the subdomains. These 
attributes facilitate the use of higher localized fidelity 
for increased solution accuracy and the ability to 
capitalize on the advantages of one mathematical 
approximation compared to another. The method is 
quite similar to the subdomain collocation approach in 
which a domain is divided into subdomains, and the 
residual is made orthogonal to a weight function and set 
to zero in an integral sense over each subdomain. Here, 
methodology is presented formulating the general 
method of weighted residuals for multiple domains by 
considering the Poisson equation, Eq. (1) for a two- 
dimensional domain for the field variable, u. 

For simplicity, the multiple-domain formulation is 
presented for two subdomains, Qj and TT (see Figure 
1(b)). Independent approximations and weight 
functions are assumed in each of the subdomains and 
continuity conditions are used to provide for a 
continuous solution across the subdomain interfaces. 
Thus, Eq. (1) is satisfied in each subdomain, 
independently, i.e., 

-V-(k i Vu i )=Q i in Q,- 


or for uniform constant, k, in each subdomain 
— k:V 2 U : = Q: ill fl • 

l l z-'l l 


subject to the boundary conditions on the subdomain 
boundaries, T,. 

As a result of the subdomain modeling, the solution 
of the problem involves an interior surface interface 

boundary, 1/ , and the information transfer across the 
boundary. Hence, the boundary surface for the z th 


subdomain is given by fj = Tf + F,' + r/ . The 
boundary conditions for the / th subdomain may be 
written as 

m, —Uj= 0 on Yf and k t — q. = Q on p. s 

d« 

The residual, R t = -k l 'V 1 u l - Q t , for each domain is 

orthogonalized by a set of weight functions, 4>, which 
yields 


j(-* i V 2 « I .-j2 / )®ldfl i =0 

£2, 

where the approximate solution, m), is sought. Using 
the general form outlined in reference 22 for the single- 
domain formulation and considering an approximate 
solution, iij, that automatically satisfies the essential 
boundary conditions, we may write the general integral 
form of the differential equation governing the potential 
flow for subdomain i as 




~kj 


f 2 ~ t 2 ~ 

O U: o u 


( 


r/ 


d;7 


d« 


dx dy 

V; 5 = 0 


-Qi 


dft ; + 


' 9i 


J 


( 2 ) 


d U: dli: dU: 

where = n , -i n v and the residual in the 

dn dx ‘ dy yi 

satisfaction of the natural boundary conditions is 
orthogonalized by a secondary set of weight 
functions <I> . 


The order of differentiation on the primary variable 
in the integral equations, Eq. (2), is reduced to obtain 
the weak formulation. In addition, for simplicity in the 
subsequent development, the tilde, (~), denoting the 
approximated values of the primary variables, u i and m 2 , 
will be omitted. Utilizing the divergence theorem, Eq. 
(2) can be rewritten, for subdomain i, yielding 


J K 


du i rHI> dli, rif|> 


dx dx dy dy 


mu- 


| kj ^-n +^-n 0 ; dr, - 

r ' dx Xi dy Vl 


jq,< i>,da + j 

rf 


k ^-a 
K i , Vi 

an 


W>dl? =0 


(3) 


Recall that the boundary integral on r„ may be 
expressed as the sum of boundary integrals associated 
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with the specified primary variable, r/' , specified 
secondary variable, F-' , and the interior surface 
interface boundary, r .with®; =0 on V/’ . Moreover, 
since®, and <I> , are arbitrary, they may be chosen, 
such that, ®; = <I> . Thus, Eq. (3) can be rewritten, for 
subdomain i, as 


J k i 

n, 


3m, 

dx 

V 


3® 3m, 

— L + — 
dx 3 y 


3®; 
3 y 


dkl i - J k i -^®;dr I 

j r i d« 


• (4) 


- \qftdr?- jG,®,dn, =o 

rf n i 


To enforce interdomain continuity on the primary 
variable, a boundary condition on the boundary, T 1 , is 
used by introducing a third approximation field in terms 
of an additional primary variable, v, (see Figure 1(b)) 
which gives rise to the continuity requirement 


v — u l = 0 on r 1 and v-n 2 = 0 on T 1 . 


These constraints can be satisfied in the integral sense 
as 


JAjG'-itJdr^O on r 1 


and 


J/l 2 (v-M 2 )dr I =0 on r 1 ( 5 ) 

r 1 

where A x and A 2 are Lagrange multipliers or weight 
functions related to the secondary variable along the 
interface. To enforce reciprocity of the secondary 
variable ( i.e fluxes or tractions) along the common 
subdomain boundary, an additional continuity 
requirement is specified. These secondary variables, 
q x and q 2 , are assumed to be independent of each 
other. These independent approximations give rise to 
continuity requirements along the interface 

q x +q 2 =0 on T 1 . 

These constraints can be satisfied in the integral 
sense as 


\A(q x + q 2 ) dr 1 =0 on T 1 (6) 
r 1 

where A is a Lagrange multiplier or weight function 
related to the primary variable along the interface. The 
integral form of Eqs. (4), (5), and (6) provides the basis 
for the subsequent spatial modeling approximations. 
The multifunctional approach developed provides a 


mechanism for coupling subdomains that have been 
discretized with diverse spatial modeling approaches. 
Of the many spatial modeling approaches, this paper 
will focus on the finite element and the finite difference 
methods and their associated coupling using the 
multifunctional approach developed. 


Eqs. (4), (5), and (6) are used to provide the 
mathematical basis for the development. In previous 
work by Aminpour et al. 10 , a similar formulation based 
on the principle of minimum potential energy is 
developed and subsequently implemented in the form 
of an element 20 . In that work as is the case in this 
study, the interdomain interface boundary is discretized 
with a mesh of evenly-spaced pseudo-nodes (open 
circles in Figure 2) that need not be coincident with any 
of the interface nodes (filled circles in the figure) of any 
of the subdomains. In this paper, the methodology has 
been extended further to encompass diverse types of 
spatial modeling methods as is indicated conceptually 
in Figure 2. 


Pseudo-nodes 



Finite element 
modeling 


Figure 2. Interface Definition. 

The generalized element equations may be 
obtained by introducing the continuity requirements 
into the weighted residual statement. Eqs. (4), (5), and 
(6) can be rewritten over an element domain as 


. , f 3 u, 3® 3m, 

\ k i — — '- + — - 
,e dx dx 3v 
■i V 


4®, 

3y 


d£^- R®,- dr 1 ' 

i r 1<? 


- R®, dr/ e - jQ^do;' = o 
rf V 

and 


JA, (v-M^dr 1 " =0, 

r 1 * 

JA 2 (v-M 2 )dr I =0, 

\X{q x +q 2 )&T lL =0. 

l e 


(7) 


Note that in the potential energy formulation 10 , the 
continuity of the secondary variables is satisfied 
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through the subsidiary conditions obtained through the 
minimization of the potential energy. In this weighted 
residual formulation, the continuity of the secondary 
variables is satisfied in a weighted residual sense and 

the Lagrange multipliers. A/ and A , are represented by 
weight functions related to the secondary and primary 
variables, respectively. 

The form of the equations for finite element and 
finite difference field approximations differs by the 
form of the element shape functions and the 
approximation selected for the weight functions, <5. For 
the generalized element expansion of subdomain /, the 
independent approximations for the element primary 
variables, (i.e., displacements or velocities), interface 
secondary variables {i.e., tractions or fluxes), the weight 
functions associated with the primary and secondary 
variables, and the interface variables are, respectively 

u, = N,u e . ; q, = R,« , ; X , = R, ; A = T 

and v = TU] 

where a ; - is a vector of unknown coefficients 


or symbolically. 


and 


+ k . s ,a, =f e . 

k Pl % +k pU I =0 

l(ki,a i )=0 


where, for 1=1,2, 


K, = I *, 

Of 


f dN J aa> + 0 n J aa> A 


mu 


dx dx dy dy 

v y 

k Pi . =- jR^Nydr 1 ' , 
r 1 " 

k . =-jafR i dr 1 *, k, = jT T R,dr r , 


and 


U, = Ka-dfl? + dr/ 


Moreover, integration over the common subdomain 
boundary, T 1 , is considered only for element edges 
along that boundary. 


associated with the secondary variables, q, , and N,, R,, 
and T are matrices of interpolation functions for the 
element primary and secondary variables, and the 
primary variables along the interface, respectively. The 
interpolation functions in the matrices R, are assumed 
to be constants for linear elements and linear functions 
for quadratic elements. The interpolation functions in 
the matrix T are cubic spline functions. Substituting 
these approximations into Eqs. (7) yields integral 
equations in terms of the weight function, <I>, which are 
given by 


Assembling the element equations over the entire 
domain, enforcing continuity of the primary and 
secondary variables only within each subdomain and 
assembling the contributions along the element edges 
on the common subdomain boundary, and noting that 
u e and u t , and f C] and f , and a! and oti are 

completely uncoupled, yields the system of equations 
given by 


( 3ivtT 


I k i 


K r , dr 1 " 

rl e 





K, 

0 

0 

K s, 

0 

V 


f. 




0 

K 2 

0 

0 

s 2 

u 2 


*2 

rKl> | <5N/ 3<I> ^ 


- 

0 

0 

0 

K h 

K, 

'2 


» = < 

0 

dx dy dy 

“ J 

e i 

K n 

0 


0 

0 

ai 


0 

i - 'T' . _ 


r e 

0 

K P2 


0 

0 

«2. 


0 


«, = mar? 


and 


jRi T N, dr 1 


jr 2 t n 2 dr' e 


lUq + 


+ 


jR, T Tdr ie 


jr 2 t t dr’ 6 



1 

0 

1 

u 


f 

or symbolically 

o 

o 

• Ul 

» = < 

0 


_K p Kj 0 

a 


0 


U] = 0 


|uj = 0 


|T T R 1 dr lC 

aj + 

|t t r 2 dr f 

\ 


_r ie 


a 2 = 0 


( 8 ) 


where K, u, and f are the assembled stiffness matrix, 
vector of primary variables and force vector for the 
entire structure, and K p , K s , K : , U|, and a are the 
assembled K„ , K s , Ki , iii, and a, for all interfaces. Eq. 

m i i 1 

(8) is obtained from the individual weighted residual 
expressions over each of the subdomains and the 
constraint integrals. 
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Multiple-Domain Modeling - Homogeneous 
Discretization 

In the context of this work, homogeneous 
discretization approaches make use of a single 


k e, = I k \ 


f W{ 3N, | 3N F dN, A 
dx dx dy dy 


mi 


and 


discretization method among all subdomains in which 
the domain is subdivided. For the finite element 

k — k 

3 2 N f 

+ 5 2 Nl 

> 

development, the weight functions are taken to be the 
finite element shape functions (;.e.,<I> = N, ). For the 

K e 2 ~ k 2 

dx 2 

V 

x=Xj dy 2 

y=yi 

"\ 

II II 

H ^ 


finite difference development, the weight functions are 
taken to be the Dirac delta function ( i. e . , 

(J> = d' [x — x, , v - v, ) ). Thus, stiffness matrices, 
k e and k s , and force vector, f e , are given in terms of 

the weight functions. The form of the coupling element 
matrices, k p and k, , that are not in terms of the 

weight functions are independent of the method of 
discretization. That is, 


k S j =- JN?Ri dr F and k Sj =-R 2 {x i ,y i ) , (9) 
r l€ 

f el = jN?adflf+ jNfodlY* 


and 


f e2 =02 (*/.?; )+?2(*,-.:y/)> 


and for the two domains, 1=1,2, 


k p =- Jr/N, dr r and k I; = |T t R,. dr r 

r 1 ' r l ‘ 

are of the same form for the finite element and finite 
difference discretizations. However, since the element 
shape functions, N„ differ for the two methods, the 
interface matrices, k p , in general, are not identical. 

Moreover, in the finite element development, at the 
element level k s = k p , and at the global system level 

K g ,=K£. By contrast, in the finite difference 
development, at the element level k s ^ k p , and at the 
global system level K s ^K p . Hence, for the finite 

difference method, the global system level matrices are 
not symmetric. The multifunctional derivation is 
general as it allows for the coupling of the primary 
variables to an independent approximation. This 
attribute is particularly important in the heterogeneous 
discretization approach described in the next section. 

Multiple-Domain Modeling - Heterogeneous 
Discretization 

Heterogeneous discretization approaches make use 
of different discretization methods for at least two of 
the subdomains in which the domain is subdivided. 
There are many combinations of spatial modeling 
approaches; however, this work focuses on the coupling 
of the finite element and finite difference methods. The 
multifunctional weighted residual formulation of Eq. 
(8) is used. Considering two domains, shown in Figure 
1(b), subdomain 1 is assumed to be discretized using 
the finite element method, and subdomain 2 is assumed 
to be discretized using the finite difference method. 
The set of element matrices becomes a hybrid of the 
matrices from the finite element method and the finite 
difference method, and 


k p =- jR F N,dr r and k I; = JT t R,. dr r . 
r 1 ' r 1 ' 

The multifunctional modeling approach has been 
generalized such that it is applicable to both 
homogeneous and heterogeneous discretization 
approaches. 

Computational Implications 

Computational implications are discussed for the 
generalized system of equations, Eqs. (8). Implications 
specific to a discretization approach are highlighted, 
where appropriate. The assembled stiffness matrix K is 
a block diagonal matrix containing the stiffness 
matrices K, of each of the subdomains along its block 
diagonal. The interface “stiffness” matrix thus contains 
coupling terms that augment the stiffness matrices of 
the subdomains along the interface. Due to the use of 
Lagrange multipliers in the constraint conditions, the 
systems are neither banded nor positive definite. 
Therefore, standard Cholesky solvers can not be used, 
unless pivoting is performed to obtain the solution. In 
addition, due to the generalization for the finite 
difference approximations, the system of equations is 
not necessarily symmetric due to different off-diagonal 
submatrices, K p and K s . The system unknowns in Eq. 
(8) consist of both primary and secondary variables 
given by the potential function, m, and the secondary 
variable coefficients, a, respectively. Generally, the 
terms in the coupling matrices, K s , are of the 

magnitude of the length of the interdomain boundary, 
which results in a marked difference in the magnitude 
of the off-diagonal terms of the system matrix 
compared to its diagonal terms. This characteristic 
produces an ill-conditioned matrix whose solution can 
cause numerical difficulties for some general-purpose 
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solvers. Hence, the coupling matrix should be scaled 
such that its terms are of the same magnitude as the 
subdomain stiffness. The upper diagonal submatrix 
blocks contain uncoupled subdomain stiffness matrices. 
The symmetry of the subdomain matrix is determined 
by the choice of the weight function, <f>. For the finite 
element discretization, the subdomain matrices are 
symmetric. However, due to the elimination of 

fictitious nodes for the imposition of boundary 
conditions and loads in the finite difference 
discretization, the subdomain stiffness matrices, K ; , 
generally are not symmetric, but they are positive 
definite and sparse. The coupling is accomplished 
through the introduction of the terms in the matrices, 
K p , K s , and Kj. The number of additional degrees 

of freedom associated with the interface is generally 
small in comparison with the total number of degrees of 
freedom in the subdomains. Thus, modeling flexibility 
is provided at a relatively small computational expense 
that may be reduced additionally as the efficiency of 
new solution algorithms for the system of equations in 
Eqs. (8) is increased. 


Numerical Results 


The multifunctional methodology for the two- 
dimensional scalar-field problem is demonstrated on a 
verification patch test case and a benchmark problem 
( i.e ., a torsion problem). The applications are 
described, and the associated results and salient features 
are discussed. Finite difference and finite element 
solutions for single- and multiple -domain 
configurations are presented to provide benchmark 
solutions for the multifunctional approach using 
homogeneous and heterogeneous spatial discretizations. 
The finite element models use four-node Lagrange 
isoparametric finite elements, and the finite difference 
models use a five-point template to approximate the 
governing differential equation. A general-purpose 
finite element code, COMET/AR 23 , is used to generate 
the finite element stiffness matrices. The mathematical 
computing program MATLAB® is used to generate the 
finite difference matrices and the interface coupling 
matrices and to solve the resulting system of equations. 

Patch Test Problems 

The fundamental concept of the patch test for the 
scalar-field problem herein is to subject a domain to 
boundary conditions that engender a linear or quadratic 
primary variable field and a constant or linear 
secondary variable field throughout the domain. For 
the governing differential equation of the form of Eq. 
(1), boundary conditions that serve this purpose are: 


i) Specified primary variable oil T which emanate 
from a linear potential field as 

u = a^x + a 2 y + a 0 
or a quadratic field as 

u = a 1 (x 2 - y 2 )+ a 2 x + a 2 y + a Q 
where a ,a ,a .and a are arbitrary constants. 

1 2 3 0 J 

ii ) Specified constant or linear secondary variable on 

r s 


q=b l x + b 2 y + b Q 

where b , and b Q are arbitrary constants. 


Consider the solution for the primary variable, 
u(x,y), in a rectangular domain (see Figure 3(a)) with 
boundary conditions of the forms indicated that yield 
the exact solution. The problem is given by Laplace’s 
equation for a planar domain as 


d^u 


= 0 , 


0 <x<a, 0 < y <b . 


a =f(x,y) 




(b) Spatial Discretization 

Figure 3. Two-Dimensional Rectangular Domain 
Geometry and Spatial Discretization. 
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Specified boundary conditions representing linear, 
bilinear, and quadratic potential fields are applied to the 


square domain. 

These boundary conditions are given as 

Linear field: 

«( 0, }') = 

2, 

u(a, y) = a + 2, and 


q„(x, 0) = 

-~<ln 

(x,b) = 0 

Bilinear field: 

h (0 ,y) = 

y> 

u(a,y) = a + y, and 


q„(x, 0) = 

=-i. 

and q n (x, b ) = 1 

Quadratic field: 

m ( 0 , y) = 

-y 

2 / x 2 2 

, u{a,y) = a —y , 


q n (x,0) = 

= 0, 

and q n (x, b) = —2b 


Single-domain and multiple -domain spatial 

modeling approaches are used to analyze the problem. 
Two single-domain models are developed each using 
either finite element or finite difference discretization 
(denoted by SD/FE and SD/FD) over the entire domain. 
Two multiple -domain models are developed using 
homogeneous modeling each using either finite element 
or finite difference discretization (denoted by MD/FE 
and MD/FD) in each of the subdomains. For 
convenience, the finite element and finite difference 
discretizations have the same number of nodes or grid 
points in each of the coordinate directions in the 
respective subdomains. A multiple-domain model 
using heterogeneous modeling (denoted by MD/HM) is 
developed with finite difference discretization used in 
one subdomain and finite element discretization used in 
the other subdomain. For the spatial modeling 
approaches, the syntax (m x n) denotes a mesh with m 
grid points in the x-direction and n grid points in the y- 
direction. For the single-domain modeling approaches, 
the entire domain is discretized using a (5 x 5) grid 
point mesh. 


For the multiple-domain analyses, a (3 x 5) grid point 
mesh is used in the less refined (coarse) subdomain and 
a (5 x 9) grid point mesh is used in the more refined 
subdomain as shown in Figure 3(b). Hence, the more 
refined subdomain has a grid point spacing in the x- and 
y-directions that is half the size of the spacing used in 
the coarse subdomain. 

Results of the analyses performed have been 
compared to the exact solution for the given boundary 
conditions and are summarized in Table 1 using 
normalized values. A value of unity implies perfect 
agreement with the reference solution. Values above 
and below unity indicate error in the computed solution. 
For boundary conditions consistent with linear and 
bilinear potential fields, the computed potential and flux 
results are exact for all analysis types. For boundary 
conditions consistent with a quadratic potential 
function, the error in the computed potential and flux is 
approximately 3% for the multiple-domain 
homogeneous finite element (MD/FE) spatial modeling, 
and the error is approximately 1% for the multiple- 
domain heterogeneous modeling (MD/HM) with 
combined finite element and finite difference 
discretization. For the given boundary conditions and 
element configuration (i.e., square or rectangular 
elements), the single-domain finite element (SD/FE) 
model reproduces the exact solution using the bilinear 
finite element. However, for a general element 
orientation (i.e., quadrilateral elements), the bilinear 
element used does not reproduce the exact solution. 
Moreover, for this problem in the multiple-domain 
analysis, error is introduced when combining finite 
element models of different discretization along the 
common boundary. However, no error is introduced 
when combining finite different models of different 
discretization along the common boundary. 


Table 1. Results of the Multifunctional Approach for the Patch Test Problems. 



Normalized Potential Function, u 

Normalized Flux, q x 

Analysis 

Order of Potential function 

Order of Potential Function 

Type* 

Linear 

Bilinear 

Quadratic 

Linear 

Bilinear 

Quadratic 

SD/FE 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

SD/FD 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

MD/FE 

1.0 

1.0 

1.03 

1.0 

1.0 

1.03 

MD/FD 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

MD/HM 

1.0 

1.0 

.99 

1.0 

1.0 

.99 


* SD/FE: 
SD/FD: 
MD/FE: 
MD/FD: 
MD/HM: 


Single-Domain with Finite Element discretization 
Single-Domain with Finite Difference discretization 
Multiple-Domain with Finite Element discretization 
Multiple-Domain with Finite Difference discretization 

Multiple-Domain with Heterogeneous Modeling (combined finite difference and 
finite element discretizations) 
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These results indicate that, for this problem, the 
finite difference solution better approximates the 
quadratic potential function. As such, the error 
obtained using the heterogeneous model (combined 
finite difference and finite element discretization) is 
smaller than the error obtained for the homogeneous 
finite element model. 

Representative Benchmark Problem 

The multifunctional methodology is demonstrated 
on a representative two-dimensional scalar-field 
application. This application is a second-order problem 
of solid mechanics ( i.e ., a torsion problem) that can be 
formulated in terms of one dependent variable (see 
reference 22 for additional problems). 


model is shown by the gray shaded region in Figure 
5(a) and in an enlarged view in Figure 5(b). 


</> = o 



The torsion of a prismatic bar with a rectangular 
cross-section is used to demonstrate the multifunctional 
capabilities of the developed methodology. The torsion 
problem reduces to the nonhomogeneous partial 
differential equation 


d 2 </> d 2 </> 


= -2G6 


in which the stress function, <j ) , must be constant along 

the boundary of the cross section, 0 is the angle of twist 
per unit length of the bar, and G is the shear modulus. 
The configuration of the bar is shown in Figure 4, and 
the analysis domain and the boundary conditions, are 
shown in Figure 5(a). 


(a) Entire Model 

«> = o 



» = o 


(b) Symmetric Model 



Figure 4. Prismatic Bar with Rectangular Cross- 
Section. 


Figure 5. Analysis Domains and Boundary Conditions 
for Prismatic Bar with Rectangular Cross-Section. 

The shear stresses in the cross-section are 
_ _ __<¥ 

44 dy ' ® dx ' 

At the ends of the bar, the first moment integrated over 
the cross-sectional area must equal the twisting 
moment. This requirement gives 

M t = 2 j <j) dxdy 


Due to the symmetries in the problem, only one 
quadrant of the rectangular cross-section needs to be 
considered. Moreover, this symmetric model is useful 
in verifying the application of mixed boundary 
conditions. That is, the application of boundary 
conditions in terms of both primary and secondary 
variables. The quadrant considered in the symmetric 


The analytical solution 24 for the stress function and 
twisting moment are given, respectively, by 

32 G6ci 2 ~ 1 / \( l ,- 1 )/ 2 r cosh(n«y/2a) 

</> = ~ L — d-l) 1 7 r cos 

n n=i,3,5,...n |_ cosh(n7E?/2aJ 

and 
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ti 2 a 


\ 
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Spatial Modeling of Prismatic Bar 

Analyses are performed for the case of b=2a ( i.e ., 
rectangular cross-section), where a and b are 
dimensions of the cross-section shown in Figure 5(b). 
Three levels of grid refinement are used for the spatial 
modeling, namely meshes of (6 x 6), (11 x 11), and (21 
x 21) grid points, each applied to the domain shown in 
Figure 5(b). A typical idealization for a (6 x 6) mesh of 
grid points is shown in Figure 6(a). Multiple-domain 
analyses with the spatial modeling of these three levels 
of grid refinement and with coincident nodes along the 
common subdomain boundary have been performed for 
comparison. For the multiple -domain spatial modeling 
with non-coincident nodes along the common 
boundary, the mesh discretization of the most refined 
domain is consistent with the discretization used in that 
same region for the single-domain analysis. The mesh 
in the less refined domain has half the “element” 
density of that used in the refined domain. The mesh of 
the entire domain is shown in Figure 6(b) and is 
referred to by the syntax (6 x 1 1 )/( 1 1 x 21). 



(a) Single-Domain (6 x 6) Mesh 



(b) Multiple -Domain (6x1 1)/(1 1 x 21) Mesh 

Figure 6. Spatial Discretization for One Quadrant of 
Prismatic Bar with Rectangular Cross-Section 


Twisting Moment for the Prismatic Bar 

Having found the values of the stress function, (j > , 
at the grid points in the solution domain by the 
respective spatial discretization approaches, the twisting 
moment may be found by repeated application of the 
trapezoidal rule for numerical integration. The 
computed twisting moment, M t , is normalized by the 
analytical solution, M t analytical , and is given in Table 2. 

A value of unity indicates perfect agreement with the 
analytical solution. Values above and below unity 
indicate error in the computed solution. Results in 
Table 2 indicate that all analyses are in good agreement 
with the analytical solution. The maximum error in any 
of the computed solutions is approximately 5% and is 
observed for the multiple-domain heterogeneous 
modeling analysis (MD/HM) using combined finite 
difference and finite element discretizations with the (6 
x 6) mesh of grid points with coincident nodes along 
the common boundary. The maximum error in this 
analysis is due to a combination of the discretization 
error and error in the multifunctional approach for 
combining the finite element and finite difference 
meshes. 

The error for the multiple-domain analyses using 
combined finite difference and finite element 
discretizations (MD/HM) with noncoincident grid 
points along the common boundary is less than 3%. 
The solution accuracy for each of the modeling 
methods increases as the mesh refinement increases. 
For the same number of nodes or grid points, the finite 
element discretization yields more accurate solutions 
than the finite difference discretization. 

The results obtained for the single-domain 
modeling ( e.g ., SD/FE and SD/FD) and the multiple- 
domain homogeneous modeling with coincident nodes 
along the subdomain boundary are identical or nearly 
identical (see the results for (6 x 6), (1 1 x 1 1) and (21 x 
21) meshes in Table 2). The results obtained for the 
multiple-domain heterogeneous modeling approach 
(MD/HM) with coincident grid points along the 
subdomain boundary are less accurate than 
corresponding results obtained using homogeneous 
modeling (MD/FE and MD/FD) but are in overall good 
agreement. With multiple-domain modeling using 
finite element (MD/FE) discretization and with non- 
coincident nodes, the twisting moment obtained is 
bounded by the twisting moments obtained using the 
less refined (11 x 11) and more refined (21 x 21) 
coincident meshes (see the results for the (6 x 1 1 )/( 1 1 x 
21) mesh in Table 2). For the multiple-domain finite 
difference (MD/FD) discretization in both domains 
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Table 2. Normalized Twisting Moment for Prismatic Bar. 


Analysis 

Type’ 

Normalized Twisting Moment, M t /{M t ) analytical 

Mesh Density 

(6x6) 

(11 x 11) 

(21 x 21) 

(6 x 1 1 )/( 1 1 x 21) 

SD/FE 

0.987 

0.994 

0.998 

- 

SD/FD 

0.974 

0.990 

0.996 

- 

MD/FE 

0.987 

0.994 

0.998 

0.996 

MD/FD 

0.975 

0.990 

0.996 

0.983 

MD/HM 

0.950 

0.974 

0.988 

0.975 


SD/FE: 

SD/FD: 

MD/FE: 

MD/FD: 

MD/HM: 


Single-Domain with Finite Element discretization 
Single-Domain with Finite Difference discretization 
Multiple-Domain with Finite Element discretization 
Multiple-Domain with Finite Difference discretization 

Multiple-Domain with Heterogeneous Modeling (combined finite difference and 
finite element discretizations) 


with non-coincident nodes, the twisting moment is 
slightly less accurate than the results obtained using the 
(11 x 11) coincident mesh, which is indicative of the 
error introduced by the finite difference interface 
constraints along the common boundary. For the 
heterogeneous modeling approach with coincident 
nodes along the interface boundary, the twisting 
moment is less accurate than the homogeneous 
approach with either finite element modeling or finite 
difference modeling. These results reveal the error 
introduced in the heterogeneous modeling approach for 
this problem due to the interface constraints. For the 
heterogeneous modeling approach with non-coincident 
nodes, the twisting moment is slightly more accurate 
than the (11 x 11) coincident mesh, which is indicative 
of the benefit gained (i.e., more accurate field 
approximation and interface constraint) by the 
combination of the finite element and finite difference 
discretizations. 


Concluding Remarks 

Multifunctional methodologies have been 
formulated for interfacing diverse domain idealizations 
including multi-fidelity modeling and multiple spatial 
modeling approaches. The methods, based on the 
method of weighted residuals, provide accurate 
compatibility of primary and secondary variables across 
domain interfaces as shown for several problems. The 
methodology has been described and demonstrated for 
the illustrative problems. These selected problems 
include a verification test case and a second-order 
problem of solid mechanics that can be formulated in 
terms of one dependent variable. The governing 
equation in each case is either the Laplace or the 
Poisson equation. In all cases considered, the results 
obtained using the multifunctional methodology were in 


overall good agreement with the reported analytical or 
reference solution. The methods have been rigorously 
developed for multiple-domain applications, and the 
robustness and accuracy has been illustrated, and the 
associated computational issues have been discussed. 
Multi-fidelity modeling approaches have been 
developed that include both homogeneous (i.e., the 
same discretization method in each domain) and 
heterogeneous (i.e., different discretization methods in 
each domain) discretization approaches. The finite 
element and finite difference methods and combinations 
thereof have been used in each of the discretization 
approaches. Results have been presented for the two- 
dimensional scalar-field multifunctional formulation 
using representative test problems in solid mechanics. 
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